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ABSTRACT 


A  frequency-domain  synthesis  method  for  multidimensional  maximum -likelihood 
filtering  of  sampled  data  obtained  from  seismic  arrays  is  presented.  This  procedure 
is  shown  to  possess  several  advantages  relative  to  the  time-domain  synthesis  technique. 
The  primary  advantage  is  that  the  frequency-domain  method  requires  approximately 
ten  times  less  computer  time  to  synthesize  the  filter  than  does  the  time -domain 
technique. 

The  details  of  a  direct  segment  method  for  the  spectral  matrix  estimation 
required  in  the  frequency-domain  approach  are  presented.  In  addition,  the  bias, 
variance,  mean  square  error,  limiting  distribution,  and  other  properties  of  the  spectral 
estimates  are  discussed.  The  details  of  a  Fortran  IV  computer  program  implementation 
of  the  frequency-domain  method  are  given.  The  experimental  results  obtained  by  pro¬ 
cessing  two  events  recorded  at  the  Large  Aperture  Seismic  Array  are  presented  as 
well  as  a  comparison  of  the  performance  of  the  frequency-domain  method  relative  to 
the  time-domain  synthesis  technique.  It  is  found  that  the  processed  noise  power  re¬ 
duction  for  the  frequency-domain  method  is  typically  about  two  out  of  a  total  of  20  db 
worse  than  that  of  the  time-domain  technique. 
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I.  INTRODUCTION 

A  considerable  reduction  of  seismic  noise  is  possible  by  employing  multi¬ 
dimensional  filtering  for  seismic  arrays.1  One  of  the  more  important  approaches  to 
seismic  array  processing  is  the  maximum -likelihood  method  which  is  equivalent  to  the 
minimum -variance  unbiased  estimator  technique1  This  multidimensional  filtering 
method  forms  a  single  output  waveform  which  serves  as  an  estimator  of  the  unknown 
signal  which  comes  from  a  fixed  direction. 

The  basic  assumption  in  our  analysis  of  multidimensional  filters  is  that  the 
th 

output,  X.  (t),  of  the  k  seismometer  may  be  written  as 


(1) 


where  S(t)  is  the  signal  waveform  which  is  assumed  to  be  the  same  in  each  seismometer 
and  N^(t)  is  the  noise  present  in  seismometer  k,  k  =  1, . . . ,  K.  In  writing  Eq.  (1)  it  is 
assumed  that  the  azimuth  and  horizontal  velocity  of  the  event,  or  signal,  have  already 
been  determined  with  sufficient  accuracy  to  allow  the  signal  waveforms  from  each 
seismometer  to  be  shifted  to  bring  them  into  time  coincidence.  In  most  applications, 
the  outputs  of  the  seismometers  are  given  in  sampled  form  in  which  case  Eq.  (1) 
becomes 


k  =  1,...,K 
m  =  0,±1,±2,. . . 


(2) 
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Only  the  sampled -data  multidimensional  filtering  problem  for  seismic  arrays  will  be 
considered. 

A  time-domain  technique  can  be  used  to  design  the  maximum -likelihood  filter. 

In  this  method  the  crosscorrelation  matrix  of  the  noise  is  measured  in  what  is  called  a 
fitting  interval,  usually  three  minutes  long,  immediately  preceding  the  event  and  using 
this  measured  matrix  the  filter  is  synthesized  by  employing  a  recursive  synthesis  pro¬ 
cedure.*  The  filter  designed  in  this  manner  is  optimum  in  the  sense  that  it  provides  a 
maximum -likelihood  estimate  of  the  signal,  provided  the  noise  is  a  multidimensional 
Gaussian  process.  It  is  also  optimum  in  the  sense  that  the  output  noise  power  is  min¬ 
imized  subject  to  the  constraint  that  the  signal  be  undistorted  by  the  filter. 

However,  in  spite  of  these  advantages,  the  time  domain  synthesis  method  is  sub¬ 
ject  to  several  severe  criticisms.  The  most  important  criticism  is  that  the  method  re¬ 
quires  a  large  amount  of  computer  time,  typically  about  30  minutes  of  7094  computing 
time  to  synthesize  a  filter  for  a  25-channel  subarray  of  seismometers.  In  addition,  this 
method  tends  to  develop  a  ,,supergainM  at  certain  frequencies,  inside  the  fitting  interval 
which  is  not  maintained  outside  the  fitting  interval  where  the  filtering  action  is  most  im¬ 
portant.  Thus,  there  is  a  considerable  effort  expended  in  this  method  in  achieving  a  large 
noise  reduction  in  the  fitting  interval  which  cannot  be  maintained  outside  the  fitting  inter¬ 
val.  Perhaps  another  way  of  saying  this  is  that  the  technique  is  too  sensitive  to  the 
assumption  that  the  noise  is  stationary.  Still  another  disadvantage  is  that  the  technique 
is  at  times  sensitive  to  the  assumption  that  the  signal  be  identical  across  the  array  of 
sensors.  This  drawback  manifests  itself  in  the  form  of  a  precursor,  when  two-sided 
or  symmetric  filters  are  used,  which  precedes  the  event.  This  can  be  troublesome  if 
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first  motion  Is  to  be  preserved  for  use  as  a  discrimination  criterion  between  natural 
seismic  events  and  nuclear  explosions. 

The  purpose  of  this  report  is  to  present  a  frequency-domain  synthesis  procedure 
for  symmetric  multidimensional  maximum -likelihood  filters  which  does  not  have  the 
disadvantages  of  the  time -domain  synthesis  method  just  mentioned.  The  theory  for  the 
method  has  been  presented  previously.  *  In  brief,  the  frequency-domain  method  is 
optimum  in  the  same  sense  as  the  time -domain  method  when  the  leqgth,  or  memory, 
of  the  filter  is  large.  The  greatest  advantage  of  the  frequency-domain  method  is  that 
it  requires  about  one  order  of  magnitude  less  computing  time  than  the  time -domain 
technique  in  the  synthesis  of  a  filter  for  a  25-channel  subarray  of  seismometers.  In 
addition,  the  frequency-domain  method  is  less  sensitive  to  the  assumption  of  noise 
stationarity  and  identical  signals  across  the  array  than  the  time -domain  technique. 
However,  the  output  noise  power  reduction  for  the  frequency-domain  method  is  typically 
about  two  out  of  a  total  of  20  db  worse  than  that  of  the  time -domain  technique.  This 
disadvantage,  however,  would  seem  to  be  offset  by  the  advantages  cited  previously. 

We  begin  our  discussion  of  the  frequency-domain  synthesis  method  by 
assuming,  for  simplicity,  that  the  noise  components  have  zero  mean  and  covariance 
matrix 


E{Njm  Nkn}  '  Vm>n>  =  • 


1  <  j,k  <  K 
v  -  m,n  s  v 


(3) 


where  E  denotes  expectation,  and  it  is  assumed  that  the  estimator,  or  filter,  is  to  use 
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2v+l  samples  extending  in  time  from  ~v  to  v*  It  will  be  assumed  that  the  noise  is 
stationary,  so  that 


TT 


pjk(m,n)  =  p.k(m-n)  =  J*  f-v(x)  6 


Jk 


—  TT 


jk 


-i(m-n)x  dx 
2tt 


(4) 


and 


fjk<x)  =  Z  PiV<m)6lmX 


m=-oe 


jk' 


j,k  —  1, . . . ,  K 


(5) 


is  the  sampled  cross  power  spectral  density  function,  x  =  mT,  T  is  the  sampling 

interval,  and  tu  is  the  frequency  in  radians /second. 

The  minimum -variance  unbiased  estimator,  or,  equivalently,  the  maximum - 

likelihood  estimate  of  S  ,  denoted  by  S  ,  can  be  written  as 

n*  J  vn 

K 


^vn  "  2  S  0^  \ 

m=-v  k=l 


m-Hi 


(6) 


where  the  filter  weights  satisfy  the  constraint 


K 

t  e. 


k=l 


km  ^mn’ 


m  =  -v, • . • , v 


(7) 


where 


6  =  1, 
mn 


=  0, 


m  =  n 
otherwise 


(8) 


and  it  is  assumed  that  a  symmetric  filter  is  to  be  used. 
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It  has  been  shown  that  for  symmetric  filters,  asymptotically  optimum  filter 


weights  are  given  by,  cf.  reference  1,  p.  16, 


1  . . “•  -'  +  2^" 


V  -  2v  'Ak«»  +  <-I>  Ak(n) 


n=l 


Re  A,  (n  — )  cos  mn  —  + 
k  v  v 


Im  A  (n  -)  sin  mn  -  ] , 
k  v  v 


k  =  1 . K 

rn  =  -V,  •  •  • » V 


(9) 


where 


S  y> 

AkW  =  “K - 

.  2  y*> 

J,k=l 


k  =  1 . K  , 


(10) 


and  {q  (x)}  is  the  inverse  of  the  spectral  matrix  {f  (x)} ,  j,k  =  1 . K.  If  v  is 

JK  JK 

A 

large ,  the  spectral  density  of  (S  -  S  )  is  given  approximately  by 

\£i  n 


K 


a(x)  =  t  Yj  qik(x>] 

j.k=l  J 


-1 


(11) 


We  note  the  following  symmetry  properties 


yx)  =  qkj(x)  ■ 

since 

V<X)  =  fkj(X)  • 

and 

yx)  =  y-x>  • 
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since 


fjk(x)  =  fik('x)  • 


In  addition,  the  constraint  equations  are  satisfied  since 


K 

E  8, 

k=l 


V-l 


km  =  ^  +  +  Z  cosmn^)] 


n=l 


/  !\  tt 

sin  (v  -  t)  m  - 

■  2v  +  —  r  v~  1 

sm  —  m  — 

2  v 


=  -  [l  +  2v-l]  =  1. 


m  =  0  , 


.In  . 

.  sin  ( —  m - mn) 

=  J_  r(.nm - 12—V - 1  j 

2V  lK'  .In  1 

sin  —  m  — 

2  v 


^  t(-l)m  -  <-l)m]  =0,  m/0, 


where  we  have  used  the  identity 


i  n 

l 


cos  ku  =  — 


x  sin  (n  +  — )  u 


k=l 


sin-u 


In  essence,  the  frequency-domain  synthesis  procedure  operates  as  follows.  The 

spectral  matrices  {f.,  (x)}  are  inverted  at  v  +  1  points  in  frequency,  namely  n  — , 

J-K  V 

n  *  0, 1 . .  to  yield  the  matrices  {q.,  (n  — )}.  The  filter  frequency  functions  A  (n  — ) 

jk  v  k  v 
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are  then  obtained  according  to  Eq.  (10).  At  this  point  the  filter  frequency  functions 


have  been  sampled  and  are  thus  specified  at^nly  v  +  1  frequencies.  In  order  t;o  obtain 

sin  -x 

the  filter  functions  for  all  -n  ^  x  ^  n,  - : -  —  type  sampling  functions  are  used  as 


follows 


sin  jx 


/  nn . 

1  X  rr  sm  v  (x - ) 

-  h  E.,v%) - y- 


n=-\rH 


.  1  .  nu . 

sin  -  (x - ) 

2  v 


k  “  ly . . .  |  1C  . 


It  is  easily  seen  that  the  Fourier  transform  of  these  frequency  filter  functions  lead  to 
the  filter  coefficients  given  in  Eq.  (9).  Thus,  the  first  step  in  the  frequencydomain 
synthesis  procedure  must  be  to  estimate  the  spectral  matrices  of  the  noise 
{f ..  (n  -)} ,  n  =  0 . .  inside  a  fitting  interval  known  to  contain  only  noise.  The 

JJC  V 

design  of  the  spectral  matrix  estimation  procedure  will  be  discussed  in  the  next 
section. 
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n. 


THE  SPECTRAL  MATRIX  ESTIMATION  PROCEDURE 


There  is  a  large  amount  of  literature  available  on  power  spectral  density 

2 

estimation  techniques.  A  book  by  Grenander  and  Rosenblatt  treats  many  of  the 
theoretical  problems  encountered  in  spectral  estimation.  Two  books  which  are  moti¬ 
vated  more  by  the  practical  considerations  in  spectral  estimation  are  due  to  Blackman 
3  4 

and  Tukey  and  Blackman.  These  contributions  deal  mainly  with  the  estimation  of  the 
power  spectrum  of  a  single  random  process.  A  discussion  of  the  problems  involved  in 
the  estimation  of  the  spectral  matrix  of  a  multidimensional  random  process  has  been 
given  by  Goodman^  and  Rosenblatt.^  Thus,  there  is  available  a  large  number  of  tech¬ 
niques  which  can  be  used  to  estimate  the  spectra  of  random  processes. 

It  is  possible  to  divide  the  spectral  estimation  procedures  into  two  broad 
categories,  namely  direct  and  indirect  methods.  In  the  direct  method  the  data  are 
transformed  immediately  into  the  frequency  domain  and  then  the  spectrum  is  measured 
using  the  transformed  data.  The  indirect  method  first  estimates  the  correlation  func¬ 
tion  of  the  data  and  transforms  this  to  obtain  an  estimate  of  the  spectrum.  This  latter 

3 

method  has  been  discussed  extensively  by  Blackman  and  Tukey.  Some  of  the  criteria 
which  have  been  used  to  judge  the  merits  of  a  spectral  estimation  procedure  are  the 
bias  and  variance  of  the  estimate,  which  will  be  discussed  in  detail  subsequently,  and 
the  amount  of  computer  time  required  to  obtain  the  estimate.  In  the  one -dimensional 
case  it  is  largely  a  matter  of  taste  as  to  which  method  one  might  use,  as  both  the  direct 
and  indirect  methods  yield  estimates  with  roughly  the  same  bias  and  variance  and  the 
computation  time  is  comparable  for  both  methods.  However,  in  the  multidimensional 
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case  it  has  been  recognized  that  much  less  computation  time  is  required  by  the  direct 

method  than  the  indirect  method,  and  both  methods  can  be  made  to  yield  estimates 

7 

with  approximately  the  same  bias  and  variance,  cf.  e.g.  Jones.  In  the  present  appli¬ 
cation  data  from  the  Large  Aperture  Seismic  Array  is  to  be  filtered  in  which  there  are 
25  sensors  in  each  of  the  21  subarrays.  Thus,  the  spectral  matrix  of  a  25 -dimensional 
random  process  is  to  be  estimated,  since  a  single  subarray  is  processed  at  a  time,  and 
in  order  to  keep  the  computation  time  within  reasonable  bounds  a  direct  method  of 
spectral  matrix  estimation  is  necessary. 

The  present  method  of  spectral  matrix  estimation  may  be  termed  a  direct 
segment  method.  The  number  of  data  points  in  each  channel  which  is  to  be  used  in  the 
estimation,  namely  L,  is  divided  into  M  segments  of  N  =  2v  +  1  data  points.  It  should 
be  noted  that  2v  +  1  is  equal  to  the  number  of  filter  points  of  Eq.  (9).  The  data  in 
each  segment  and  each  channel  are  transformed  into  the  frequency  domain  and  these 
transforms  are  used  to  obtain  an  estimate  of  the  cross -spectra  in  the  segment.  The 
stability  of  the  estimate  is  then  increased  by  averaging  over  the  M  segments.  We  now 
describe  the  method  in  some  detail. 

The  transform  of  the  noise  data  in  the  n^1  segment,  j^1  channel,  and  at  the  <tt^1 
frequency  is 


6imt(2TT/N-l) 


(12) 


j  =  l 
-1  =  0 


2 

M, 


K, 

N-l 


n  =  1 


9***9 
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where  w  ,  m  =  1,. . .  ,N  are  the  coefficients  of  the  weighting  function.  Let 


X(l)  =  l 


2tt 

N-  1  * 


,  =  0  Nil 

•  •  •  9  2  » 


(13) 


so  that 


-1/2 


N 


S.  (X)  *  (N)  2  w  N  u 

jn  m  j,m+(n- 


_imX 


1)N 


(14) 


As  an  estimate  for  f.,  (X)  we  take 

jk 


M 


yx>  -  b  \  y»  sl<x>  ■  - 1 . k 

J  n=l 


(15) 


The  mean  value  of  this  estimate  is 


E{fjk(X» 


1  M 

*  E  M  2,  SJn<X>Sta<X> 
n=l  J 

L  M  N 

=  E  Z-JT7  /  / .  w  w  , 

MN  u.  LJ m  m 

n=l  m,m-l 


^j,m+(n-l)N  ^k,  m'+(n-l)N  ^ 


N 


i(m-m')X 


=  ^  Y.  w  w  ,  p„  (m-m')  € 

fsg  »  m  m  'llr 


ci(m-m’)X 


m 


■J  m  m  jk 

,  m  =1 


TT 


=  s  y*>  i  wN<x-i>  i 2  %  • 

-TT  J 


(16) 


where  the  frequency  window  W^(x)  is  defined  as 
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wN(x) 


=  (N) 


-1/2 


N 


m=l 


w  € 
m 


-lmx 


(17) 


and  thus 


TT 


W 


m 


-  N1'2  J  WN(x)€ta  &  . 

-TT 


(18) 


As  usual,  we  require  that  W^(x)  be  a  reasonably  good  window  in  the  sense  that  |  W^(x)| 
approaches  a  delta  function  in  such  a  manner  that 


TT 


I  Iwm«I2  i  ■ 


—  TT 


(19) 


In  this  case  f^(X)  is  a  reasonably  good  estimate  for  f^(\). 

Since  f.,  (X)  is,  in  general,  complex  it  is  convenient  to  define 


fjk(>0  -  Cjk(X)  +  iqjk(X) 


j,k  —  1, . . . ,  K  , 


(20) 


where  cjk(X)»  qjk(X)  are  real -valued  functions  and  are  known  as  the  cospectrum  and 
quadrature  spectrum,  respectively.  Similarly, 


fJkM  -  c.kft)  +  iqjk(X)  . 


(21) 


It  thus  follows  that  c.,  (X),  q.,  (X),  are  reasonable  estimates  for  c.,  (X),  q  (XX  respectively 

jk  jk  jk  jk 

A 

The  mean-square  value  of  c^(X)  is 
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E{cjk(X)}" 


E  l 


N 


/,  A  w  w  ,  w  ,,  w 

M2N  V  .  V7  , t  m  m  m  m 

n, n-1  m,m  ,m  , 


m’”=l 


N.  ,/  i\MNt  i\mN-  . »  i\mNu  t.tj/  .  1VMcos  (m-m*  )Xcos  (m'  ,-m,n 

j,m+(n-l)N  k,m  +(n-l)N  j,m  +(n  -1)N  k,m  +(n  -1)N 


(22) 


In  order  to  proceed  with  the  analysis,  we  must  at  this  point  introduce  the  assumption 
that  {Nj,m}  is  a  multidimensional  Gaussian  process,  so  that 

E-f  N  N  N  N  \  = 

1  j,m^n-l)N  k,mNKn-l)N  j,m’NKnM)N  k,m",-Kn’-l)N; 

pjk(m-m’)  p.k(m"-m,")  +  p.j[in-m"  +  (n-n')  N]  • 
pkk*m’"m"  *  +  (n_n’>  +  Pjk(rn'm'M  +  (n'n’)  N 1  * 

p^lm’-m”  +  (n-n')  N  ]  .  (23) 


Using  the  above  result  in  Eq.  (22),  we  obtain,  after  some  manipulations, 

VAR  {£..(>.)  }  =  f  Re  /  /  {f  (x)  f  <x')  + 

J  “*  TT  “*  TT  JJ 

sin^N(x-x') 

fk(x)£  (x')}  (— — — )  [|WN(x-X)WN(x'+X)|  + 

M  sin  -  (x-x  ) 

WN(x-X)W*(x+X)WN(x’+X)W*(x'-\)]  ^  .  (24) 


)X. 
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At  this  point  it  is  necessary  to  place  some  regularity  conditions  on  W^(x)  which 

2 

guarantee  that  |W^(x)|  behave  like  a  delta  function  at  x  =  0,  as  N-“,  i.e. ,  in  addition 
to  Eq.  (19)  we  require, 

,  L  |WN(X)|2  I  - 1  ■  (25> 

lxl<6 


for  any  6,  as  N  -♦  ®  , 


|WN(x)|  -  0  , 


(26) 


uniformly  when  |x|  b  €,  as  N  -*  °°,  and 


MAX 


I  |wN<x)|2  |wN(rtxo)|2 f 

JT _ 

/  lWN<X>!4  S 


“  1 


-  o  , 


(27) 


as  N  for  any  A  >  0.  Under  these  conditions  it  may  be  shown,  in  a  manner  similar 

2 

to  that  given  by  Grenander  and  Rosenblatt  ,  pp.  137-145,  that  Eq.  (24)  becomes,  as  M, 
N 


VAR(VX»  ?^T_({f)j<x)W;l)  +  c)kW_1k<x>}  ' 

A  y 

|WN(X-X)|  ,  X^O.n 

{f]i(x)tkk(x)+ik(x)-^k<x)}  • 

|WN(x-X)|4  ,  X  =  0,  TT.  (28) 
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Similarly, 


VAR  {qjk(X)} 


1 

2MN 


;  {fjjWf^w+^w-^x)}  • 


A  fly 

!WN<x'X)1  2tt  *  X  *  °»  n 


=  0,  \  -  0,  n  , 


(29) 


since  q.R(0)  =  qjk(0)  =  0,  and  q.^n)  =  qjk(n)  =  0.  If  f..(x),  f^x),  cjk(x),  qjk(x),  are 
reasonably  constant  in  the  vicinity  of  x  =  X,  then  (28)  and  (29)  can  be  rewritten  as 


VARV)}  "SN  + 

/  |wN(x-x)|4  ,  x^o,  T, 

“  TT 

TT  4  , 

J*  |WN(x-X)|4  ^  ,  X  =  0,  TT  (30) 

“TT 

VAR  {V<X>}  “  2>5n  V’  fkk(X>  +  ^k<X>  -  CJlc<X»  • 

/  |WN<x-X)|4  fjj  ,  X  /  0,  TT 

—  n 

=  0,  \  =  0,  TT  .  (31) 
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In  the  important  case  of  uniform  weighting,  =  1,  m  =  and 


N 


|wN«r  =  s  I  E  e 

m=l 


-imx,  2 


1 

N 


.  N 
SM>  2  * 

sin  |x 


(32) 


It  is  easily  seen  that  this  weighting  function  satisfies  the  conditions  in  Eqs.  (19),  (25), 
(26),  (27).  We  note  that 


iwN«r  = 


,  N 
=  -  7 

N  L  ,  . 
m,m  =1 


-i(m-m')x 


N 


-  E  »-■>=*->  e 

m=-N 


•imx 


(33) 


so  that 


nv>i4 1=  E<i-^> 

-tt  m=-N 


=  2N  +  1  -  ^  ^  (£)  (2N3  +3N+1) 


21^  +  1  _  2N 
3N  "3 


(34) 


Thus,  using  (34)  in  (30),  (31), 
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X  ^  0,  IT 


vAR<yx»  "3b  t^w^w+^w-^w} 

4¥)fi*w+c]kw^wt '  x=0>n  (35) 

VAR  {qjk(X)>  =  3^  {f..(X>  fkk(X)  +  q=k(X)  - c^X)}  X  +  0,  tt 

=  0,  X  =  0,  rr  •  (36) 

It  follows  from  (35),  (36)  that  c  (X)  and  q  (X)  are  consistent  estimates  for  c  (X),  q  (X), 

jk  jk  jk  jk 

respectively,  and  that  the  variances  of  these  estimators  go  to  zero  as  M  -*  ®. 
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ni.  BIAS  AND  MEAN  SQUARE  ERROR  OF  SPECTRAL  ESTIMATES 


The  bias  of  the  estimate  of  the  cospectrum  is  defined  as 


bNjk(X)  E  ^Cjk(X)^  Cjk(X)  ’ 


(37) 


and  provides  an  indication  of  the  error  between  the  mean  value  of  the  estimate  and  the 
true  value,  at  each  frequency.  We  may  evaluate  the  bias  for  the  uniform  weighting 
function  by  using  Eqs.  (16)  and  (21) 


bNik^  =  Fj  2  Pik(m-m’)  cos  (m-m’)X  -  J  p.k(m)  cos  mX 

J  m,m'=l  J  m=-«  J 

N-!  I  . 

=  2  (*  -  ^  )  pik(m> cos  mx  -  2  pik(m> cos  mX 

m=-N+l  •*  m=-<=  J 


N-l 

--  E 


m=-N+l 


pik 


(m)  cos  mX 


Pjk 


(m)  cos  mX 


1  V 

'  S„So 


TT 

cos  mX  J 

-TT 


c!  (x)  sin  mx  ^ 
jk'  2tt 


,  ,  .  dx 

cjk(x)  s 


A  _  gim(x-X) 
x  _  gi(x-X) 


im(x+X)  . 

— - }  c’  (x)  — 

1  _  gi(x+X)  i  Cjk'  *  2tt 
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+ 


o  11 

S  J  < 

—  TT 


sin  (x-X)  —  sin  m(x-X)  4-  sin  (m-1)  (x-X) 

4  [  1  -  cos  (x-X)] 


sin  (x+X)  —  sin  m(x+X)  +  sin  (m-1)  (x+X)  i  ,  .  .  dx 
4  [  1  —  cos  (x+X)]  *  Cjk'X'  2tt 


tt 


^  -  p  r 

N  j  1- 


sin  (x-X) 


,  .  .  dx 


-  TT 


cos  (x-X)  jk'  '  2tt 


(38) 


where  P  denotes  that  the  principal  part  of  the  last  integral  is  to  be  taken.  Similarly, 


bNjk(x)  - E  { Vx)} '  Vx) 


f"  -2 

N  J  1- 


sin  (x-X) 


,  .  .  dx 

q,u(x)  ^  * 


—  TT 


cos  (x-X)  Hjk'  '  2tt 


(39) 


Thus,  the  bias  of  the  estimators  of  the  cospectrum  and  quadrature  spectrum  decrease 
to  zero  at  the  rate  of  N 

The  mean  square  error  of  the  cospectral  estimate  is  defined  as 


E{cjk(X)-cjk(X)r  =  b^jk(X)+VAR{cjk(X)}  , 


(40) 


which  may  be  evaluated  for  the  uniform  weighting  function  by  using  (35),  (38) 


+  +,‘ pj"  c'k(x)  f]2.  M0.m(4l) 
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We  assume  that  M  is  of  the  form  M  =  cLa,  where  c  is  a  constant,  and  it  will  be  recalled 
that  L  is  the  total  number  of  data  points  used  in  the  spectral  matrix  estimation  and  is 
equal  to  MN.  We  will  determine  the  c  and  a  which  ensure  that  the  mean  square  error 
is  of  smallest  order  as  M,  N  -*  ».  it  is  easy  to  see  from  (41)  that  a  =  2/3  and 


c 


_  r  VX)  +  ^ ,1/3 

TT  •  ,  .v  ,  o 


f= 


sin  (x-X) 


—  TT 


,  .  v  dx.2 

cos  (x-X)  Cjk'X'  2tt 


X  /  0,  TT  , 


=  r  fjj(x)  WX)  *  °jk(X)  qjk(X)  Y !/3  ,  X  =  0,  TT 

1  TT  .  /  *  v  J  r 


sin  (x-X) 


—  TT 


,  n  dx  2 
cos  (x-X)  CjkW  2tt 


(42) 


Similarly,  for  the  quadrature  spectral  estimate  a  =  2/3  and 


=  { 


il/s 


TT 


«ipf 


“TT 


sin  (x~X)  Q.  (x)  dx  2 
cos  (x-X)  qjk'X'  2tt  J 


X  /  0,  TT 


(43) 


-2/3  2/3 

Thus,  the  mean  square  error  approaches  zero  at  the  rate  L  '  and  M  =  cL  , 
_  1  1/3 

N  =  c  L  .  The  above  formulas  are  not  too  useful  in  designing  the  spectral  matrix 
estimation  procedure  since  the  spectra  are,  of  course,  unknown.  In  addition,  as  far  as 
designing  a  multidimensional  maximum -likelihood  filter  is  concerned,  the  mean  square 
error  may,  or  may  not,  be  the  important  criterion.  The  merits  of  a  spectral  matrix 
estimation  procedure  must  be  determined  by  using  digital  computer  experimentation. 
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However,  the  above  results  are  useful  in  the  sense  that  they  indicate  that,  as  L  is 
increased,  M  should  be  increased  at  a  faster  rate  than  N  if  the  mean  square  error  of 
the  estimate  is  to  be  minimized. 
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IV.  LIMITING  DISTRIBUTION  OF  SPECTRAL  ESTIMATES 


If  the  regularity  conditions  for  W^(x)  given  in  (19),  (25),  (26),  (27)  are  satisfied, 
and  if  we  assume  that  {N^}  is  a  multidimensional  Gaussian  process  in  order  to  avoid 
some  complicated  regularity  conditions  which  make  the  following  result  valid  for  non- 
Gaussian  processes,  then  a  straightforward  extension  of  a  limit  theorem  of  Rosenblatt^ 


shows  that  the  c.,  (X),  q.,(X),  j,k  =  1, . . .  ,K  are  jointly  asymptotically  normally  distribut< 
jk  jK 

with  variances  and  covariances  given  by 


TT 


limit  (2MN)  {  J  |WN(x-X)|4fi  J  |Wn<*-h)|*  } 


TT 


4  dx  i-l/2 


M,N-« 


-tt 


—  tt 


COV  {  SaSW,  £y6<U)}  -  Re  [f*YW  +  f*6W  fBV(X)] ,  X  *  0,n 


=  2  V  WX)  WX)  +  fa6<X>  V”  1  •  X  =  •  <44> 

limit  <2MN){  J  |WN(x-X)|4^/  |WN<x-u>|4  §|F1/2COV  {c  (X),  q  (X)} 

M,N-*»  —  n  —  tt 

=  6.  Im  [f*  (l)f  .(\)-f*Wf0  (X)]  ,  X^O,  tt 
Xp.  ay  36'  a6  3y 

=  0  ,  X  =  0,  tt  (45) 
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limit  (2MN)  {  J  |WN(x-X)|4g  f  |WN(x-U)|4  %}'1/2  COV  {5  IX),  q  (X)} 

M,N-*®  —  rr  —  tt 

=  V  Re  I  'aY<l)  f -  f<*6W  f ?v(x)  1  •  M  °>  " 

=  0  ,  \  =  0,  TT  .  (46) 

It  is  possible  to  use  these  formulas  to  establish  asymptotic  confidence  intervals  for 
the  spectral  estimates.  A  different  approximation  for  the  asymptotic  distribution  of  the 
spectral  estimates  has  been  given  by  Goodman^  in  terms  of  the  complex  Wishart 
distribution. 
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V.  NONNEGATIVE -DEFINITENESS  AND  INVERSION  OF  ESTIMATED  SPECTRAL 
MATRICES 

We  obtain  from  Eqs.  (14)  and  (15)  that 


K  a  M  K  N 

Y.  a.a*  f..  (x)  =  —  /  /.  T.  w  a.  N.  .... 

.  ]  k  ik  L  u.  .u.  u  .  mi  i,m+(n-l)N 

j,k=l  J  J  n=l  j=l  m=l  j  j»  \  / 


,imx 


—  0,  “TT  —  X  —  TT,  (47) 


so  that  the  spectral  matrix  is  nonnegative -definite  at  all  frequencies.  This  is  a  highly 
desirable  property  since  inverses  of  spectral  matrices  must  be  computed  in  the 
frequency-domain  synthesis  procedure.  If  the  spectral  matrix  is  singular  at  a  particular 
frequency  then,  physically  speaking,  this  means  that  infinite  suppression  of  the  noise  at 
that  frequency  is  possible.  As  a  practical  matter,  it  is  not  possible  to  achieve  infinite 
suppression,  so  that  the  spectral  matrix  will  not  be  singular.  This  conclusion  has  been 
borne  out  by  computer  experimentation,  as  no  singular  spectral  matrices  have  ever 
been  encountered. 

In  general,  the  estimated  spectral  matrix  will  be  nonnegative -definite  and 
Hermitian.  Let  us  denote  such  a  matrix  by  C  =  A  +  iB,  where  A,  B  are  real  matrices. 

It  is  easily  seen  that  A  is  symmetric,  B  is  skew  symmetric  and  that  neither  A  or  B 
need  be  nonnegative -definite  if  C  is  nonnegative -definite.  Let  the  inverse  of  C  be 
denoted  by  D  +  iE.  It  is  easily  seen  that 


1 

D  ,  E 

i 

A  1 

B 

r 

-  E  1  D 

1  J 

L‘Bi 

A 
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Thus,  the  inverse  of  the  matrix  C  may  be  obtained  by  inverting  the  matrix 
The  advantage  of  performing  the  matrix  inversion  in  this  manner  is  that  this  matrix 

g 

is  invertible  if  and  only  if  the  matrix  C  is  invertible,  cf.  Goodman.  It  is  for  this  reason 
that  this  method  of  inverting  the  estimated  spectral  matrices  is  used  in  the  computer 
program  implementation  of  the  frequency-domain  synthesis  procedure. 
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VI.  COMPARISON  OF  DIRECT  SEGMENT  METHOD  WITH  OTHER  METHODS  OF 
SPECTRAL  MATRIX  ESTIMATION 

The  amount  of  computer  time  required  to  estimate  the  crosscorrelation  matrix 
for  N  lags  is  given  approximately  by 

T  =  L  K2  N  (n  +  a)  seconds  ,  (49) 

0 

where  (a,  a  are  the  multiply  and  addition  times,  respectively,  in  seconds.  It  is  seen 
from  Section  II  that  the  amount  of  computer  time  required  to  estimate  the  spectral 
matrix  at  (N  +  l)/2  frequencies  is  approximately 


and 


T  =  L  K  (K  +  N)  (n  +  a)  seconds  , 
s 


T 

_£  KN 
Ts  ~  K  +  N 


(50) 


(51) 


If  K  and  N  are  large,  the  ratio  Tc/Tg  can  be  very  large.  If  K  =  25,  N  =  21,  then 

Tc/Ts  =  11.4.  Thus,  we  see  that  the  amount  of  computer  time  required  to  estimate 

the  spectral  matrix  is  inherently  much  smaller  than  that  required  to  estimate  the 

crosscorrelation  matrix  when  K  and  N  are  large.  An  even  greater  saving  is  possible 

9 

by  using  the  Cooley-Tukey  algorithm  to  transform  the  data  into  the  frequency  domain. 
In  this  case,  the  amount  of  computer  time  becomes 


Tg  =  L  K  (K  +  F^)  (^i  +  a)  seconds, 


(52) 
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and 


K  +  F. 


KN 


N 


(53) 


where  F^  is  the  sum  of  the  factors  of  N  —  1.  If  N  —  1  =  20  =  5  •  2  •  2,  then 
FN  =  2  +  2  +  5  =  9  and  Tc /Yg  =  15. 5. 

If  an  indirect  spectral  matrix  estimation  procedure  is  used  in  which  the  estimates 
of  the  crosscorrelation  function  are  transformed  into  the  frequency  domain,  the  amount 
of  computer  time  becomes 


Tj  =  (LK2N  +  K2N2)(n  +  a) 


—  LK^N(n  +  a)  seconds 


(54) 


=  T  ,  L  »  N  . 
c’ 


It  has  been  assumed  that  the  Cooley-Tukey  algorithm  is  not  employed  since  the  saving 

would  be  incurred  in  the  second  term  above  which  is  already  negligible  compared  to 

the  first  term.  Thus,  Tj/Tg  is  quite  large  and  Tj/T'g  is  even  larger.  This  disadvantage 

makes  the  indirect  method  undesirable  for  spectral  matrix  estimation. 

The  amount  of  computer  time  required  to  synthesize  the  maximum -likelihood 

2  3 

filter  in  the  time  domain  is  2.  5  N  K  (|j  +  a)  seconds  and  in  the  frequency  domain  is 
3  3 

2N  K  ,  where  it  is  assumed  that  8  K  operations  are  required  to  invert  a  matrix  with 
complex  elements,  according  to  the  procedure  given  in  Section  V,  cf.  reference  1, 
pp.  18-19.  Thus,  the  total  computing  time  required  in  the  time  domain  is 
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(55) 


Tt  =  (L  K2  N  +  2. 5  N2  K3)  (|j  +  a)  seconds 


and  in  the  frequency  domain  is 


Tp  =  (L  K  (K  +  N)  +  2  N  K3)  (|i  +  a)  , 


(56) 


assuming  the  Cooley-Tukey  algorithm  is  not  used.  Thus 


L  K  N  +  2.  5  N2  K2 
L  (K  +  N)  +  2  N  K2  * 


(57) 


which,  if  L  is  1800,  corresponding  to  3  minutes  of  noise,  20  data  points  per  second  and 
every  other  data  point  skipped,  and  K  =  25,  N  =  21,  becomes  T^,/Tp  =  15.  In  practice 
this  large  saving  is  not  achieved,  of  course,  due  to  indexing  operations  and  tape  reading 
which  are  common  to  both  methods,  but  savings  of  approximately  a  factor  of  ten  have 
been  realized. 
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Vn.  COMPUTATION  OF  INPUT  AND  OUTPUT  NOISE  POWER  INSIDE  THE  FITTING 
INTERVAL 


The  input  noise  power  at  the  •£t*1  frequency,  inside  the  fitting  interval,  is  obtained 
from  the  spectral  estimation  program  as 


NS 


PINW  (N-l)K  xjj 


2  , 


(58) 


and  the  total  noise  power  is 

N-3 
2 


IN 


fj—  1 


IN 


IN' 


IN  '  2 


(59) 


The  justification  for  this  definition  of  is  apparent  and  the  justification  for  defining 

Pj^  in  Eq.  (59)  is  that  P^  may  be  written  as,  after  some  manipulation, 


IN 


i  K  L  0 

=  JL  v  v  N2 

LK  H  iNjm  * 


(60) 


j=l  m=l 


which  is  numerically  equal  to  the  average  input  noise  power. 

Let  us  consider  the  power  measurement  which  is  made  by  summing  the  squares 
of  the  processed  noise  samples  at  the  midpoints  of  the  segments  used  in  the  spectral 


matrix  measurement,  i.  e. , 

N-I 

M  NS  2 


POUT  M  S  (  Z  Em.»  ®jm  Nj,m+(n-l)N^  * 
n=l  j=l  m=-(-^— )  J  J 


(61) 


We  may  show,  after  some  manipulation,  that 
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(62) 


N^l 

2 


POUT  N-l  N-l  ^  'Ijk^  1  }  * 


where  {q  (£)}  is  the  inverse  of  the  estimated  spectral  matrix  {f  ({,)}.  The  frequency- 
Jk  Jk 

domain  computer  program  uses  the  inverse  spectral  matrices  to  compute  as 

indicated  in  Eq.  (62).  It  should  be  noted  that  this  estimate  is  based  on  M  out  of  a 
possible  L  data  samples  in  the  fitting  interval  and  must  be  regarded  as  an  approximation 
for  the  true  output  power  in  the  fitting  interval  which  is  equal  to  the  sum  of  the  squares 
of  the  L  data  points  in  the  filtered  output  trace.  In  a  similar  manner,  the  output  power 
at  the  f,**1  frequency  is  computed  as 


N  .  NS 

pout^  =  n7!  FT!  2 


(63) 


j,k=l 


which  can  be  shown  to  be  equal  to  the  following 


1  M 

M(N-1)^  S 
n=l 


S  -  AkW  Nk,m+(n-l)N  € 
k=l  m=l 


-im^(2TT/N-l) 


(64) 


This  latter  quantity  is  what  is  obtained  by  transforming  the  noise  data  in  the  segments 
used  in  the  spectral  matrix  estimation,  applying  the  filter  A^(-t),  and  then  summing  the 
squares  of  the  magnitudes  of  results  in  each  segment.  This,  then,  is  the  justification 
for  the  definition  of  ^OUT^*  ^  addition,  we  have 
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N-l 


p 

OUT 


POUT^ 


(65) 


which  is  a  desirable  property  for  power  estimates.  It  is  interesting  to  compare 
Eqs.  (62),  (63)  with  (11),  as  this  provides  further  justification  for  the  definitions  used. 
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Vm.  DESCRIPTION  OF  COMPUTER  PROGRAM 


The  Fortran  IV  computer  program  implements  the  synthesis  of  the  maximum  - 
likelihood  filters  in  the  frequency  domain  according  to  the  algorithm  given  previously. 
The  essential  outputs  of  this  program  are  the  above  filters  and  SC -4020  plots  of  the 
processed  traces.  The  program  works  with  data  duplicated  from  Large  Aperture 
Seismic  Array  tapes  on  a  FORTRAN  compatible  tape  known  as  EDIT-3  tape.  A 
description  of  the  program  is  as  follows. 

INPUT  TO  PROGRAM 
Card  1  Format  (815) 

NSKP,  WARMUP,  BLKSZ,  NSAMB,  NIM,  NS,  NFP,  POCTR 
where 

NSKP  is  the  number  of  EDIT -3  data  records  or  number  of  data  samples  skipped  initially. 
WARMUP  is  the  number  of  EDIT-3  records  used  to  determine  average  D.  C.  levels  for 
each  seismometer.  The  block  over  which  the  average  D.C.  is  estimated 
immediately  preceeds  the  fitting  interval. 

BLKSZ  is  the  number  of  points  of  data  used  in  each  sample  block. 

NSAMB  is  the  number  of  sample  blocks. 

NIM:  Every  NIM^1  point  is  used  as  data,  in  the  fitting  interval,  the  others  are  skipped. 

It  follows  that  the  total  number  of  EDIT-3  records  spanning  the  fitting  interval  is 
BLKSZ  *  NSAMB  *  NIM  which  amounts  to  (BLKSZ  *  NSAMB  *  NIM  -  l)/20  seconds  of 
sampled  data. 

NS  is  the  number  of  seismometers. 

NFP  is  the  number  of  filter  points  and  is  an  odd  integer.  At  present  we  make  NFP  and 
BLKSZ  equal. 

POCTR  is  a  printout  control  parameter.  If  it  is  not  zero,  the  cross  power  spectral 

density  matrices  will  be  printed  out.  If  it  is  zero  or  left  blank,  these  matrices 
will  not  be  printed  out. 

Cards  2  to  2  +  NS- 1  Format  (15,  5X,  A6,  4X,  2F10. 3)  contain  HCH(i),  ISEIS(i),  XC(i), 
YC(i),  i  =  1,  NS  where 
HCH  is  the  LASA  seismometer  number. 

ISEIS  is  the  alphabetic  name  of  the  seismometer  (such  as  A0,  B3). 

XC(i)  is  the  N.S.  subarray  projection. 
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YC(i)  is  the  E.W.  subarray  projection  of  the  seismometer  numbered  IICH(i). 

Card  2  +  NS  Format  (2F 10. 3) 

AZ  =  azimuth 

HVEL  =  horizontal  phase  velocity 
Card  3  +  NS  Format  (215) 

NJUMP  =  number  of  EDIT-3  data  records  skipped  before  plotting.  NSQSM  =  No.  of 
records  processed. 

Card  4  +  NS  Format  (F10. 3,  15) 

TRCMX.  The  positive  scaling  for  plot. 

Total  scaling  range  is  2  TRCMX. 

NSITE  LASA  Site  Number. 

The  input  tape  is  an  EDIT-3  tape  of  the  event.  The  following  constraint  must  be  obeyed 
by  the  data: 

BLKSZ  *  NIM  <  66 
NFP  *  NIM  <  66 

and  BLKSZ  (and  NFP)  must  be  odd. 


OUTPUT 


SC-4020  hard  copy  plots  of  the  center  seismometer,  WDS,  DS,  FS,  and  measured 
variances  of  FS,  WDS  and  average  power  of  all  seismometers  taken  over  blocks  of  100 
consecutive  data  points.  FS  denotes  filter  and  sum  and  is  the  designation  for  the 
maximum -likelihood  filter  which  employs  NFP  filter  points,  WDS  denotes  weighted 
delay  and  sum  and  designates  a  one-point  maximum -likelihood  filter,  i.  e. ,  NFP  =  1, 
and  DS  denotes  the  sum  of  the  delayed  data. 

TAPES 


EDIT-3  hang  on  B8 
Buffer  on  A7. 

However,  this  may  be  changed  at  load  times  by  changing  $  ATTACH  and  $  AS  card 
contents. 

RUNNING  TIME  ON  7094 

For  WARMUP  =  600,  NS  A  MB  =  90,  NIM  =  2,  BLKSZ  »  NFP  *  21,  NSQSM  =  7200, 
6  minutes  of  data,  the  7094  computer  running  time  is  about  10  minutes,  consisting  of 
4  minutes  to  synthesize  the  filter  and  6  minutes  to  obtain  all  the  processed  traces. 
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JOB  STACKING 


If  several  sets  of  data  cards  are  run  with  one  program  deck,  then  the  jobs  will 
run  consecutively  without  intervention  and  there  is  but  one  setup  for  tapes. 

A  listing  of  the  program  is  given  in  Appendix  A. 

A  detailed  description  of  the  flow  of  the  program  is  as  follows.  Note  that  the 
uniform  weighting  function  was  used  and  that  the  Cooley-Tukey  algorithm  was  not 
employed.  However,  a  later  version  of  the  program  will  use  this  algorithm. 

1.  Compute  the  delays  for  the  NS  channels  corresponding  to  the  azimuth  and 
horizontal  velocity  of  the  event. 

2.  Apply  the  computed  delays  to  the  data,  and  let  the  delayed  data  sample  in  the  j1*1 
channel  at  time  point  m  be  denoted  by  X^m,  to  be  consistent  with  previous  notation. 
Let  m  =  1  designate  the  beginning  of  the  warmup  period. 


3. 


j  =  1, . . .  ,NS,  where  L’  =  WARMUP,  which 


usually  consists  of  600  data  points  immediately  preceding  the  fitting  interval. 
No  data  samples  are  skipped  in  this  computation. 


4.  Let  X!  =  X.  —  X.  and  compute 


jm  jm  j 


€im^(2n/NFP-l) 
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7. 


Compute 


NS 

£  V0 

■  -£ — 


y  Kf  Oj  •  •  •  | 


NFP-1 


Yj  Qr.j('t')  —  !»•••»  NS. 

j,k=l  J 


8. 


Note  that,  for  each  l,  the  denominator  need  be  computed  only  once  for  all  k. 
Compute 


1  m  NFP-1 

0,_  =  [AJ)  +  (-l)mAt(^)  + 


jm  NFP-1  k 
NFP-3 


k'  2 


2  Yj  Re{Ak(^)}  cos  (nrt,  ^pp_1  )  + 
<t=l 


Im  {Ak(^sin  (m-l  j^|p7y 


)]  ,  j  =  1 . NS 


„  ,  NFP-1, 

m  =  0,  ±1, . . .  ,±  (  — o — )  . 


9.  Compute  the  zero -lag  correlation  coefficients 

NFP-1 

2 

p  (0)  =  (NFP-1)'1  Y  f  ,  j  =  1, ....  NS 

3  ,  .NFP-3,  JJC 

b=~( — 2 — )  k  -  J  • 

Set  p  (0)  =  p  (0),  k  <  j,  j  =  1, . . .  ,NS.  Note  that  p  (0)  is  numerically  equal 
jk  kj  jk 

to  the  estimate  obtained  by  correlating  data  directly,  namely 
1  L 

^jk(0)  =  L  m^1  Xj,m(NIM)  +  L’  Xk,  m(NIM)  +  L*  * 


10.  The  filter  coefficients  for  WDS,  0k>  k  =  1, . . .  ,NS,  are  obtained  by  inverting 
the  zero -lag  correlation  matrix  bordered  by  zeroes  and  ones,  just  as  in  time- 
domain  method.  A  Lagrangian  multiplier  is  also  obtained  in  this  computation 
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which  is  numerically  the  same  as  the  output  power  of  WDS  in  fitting  interval, 
PWDS* 


11.  Compute 


1  NS  . 

PIN^  =  NFP  —  1  NS  2  ’ 

J=1  JJ 


t  =  0,.. 


NFP-1 
*  2 


and 


NFP-3 


NFP  —  1 

P  =  P  (0)  +  P  (  — - 

IN  INV  '  IN  V  2 


)  +  2  Yj  PjnW  • 


12.  Compute 


NS 


POUT^ 


(NFP  -  1)^  1  2 


NFP-1 


and 


NFP-3 


NFP-1  v> 

=  %rr<0)  +  pr„rr(-T^)+  2  Z  pnm-<«  • 


OUT  OUT'  '  OUT'  2 


OUT' 


13.  Compute  signal -to -noise  ratio  gain  of  FS  at  t1*1  frequency 


Pin  W 

gfs<«  ■  10  ‘°SlO  p^lj 


and  total  gain 


gfs  =  lolQSop 


IN 


OUT 


14.  Compute  signal -to -noise  ratio  gain  of  WDS 


gwds  =  10  lo%)P 


rIN 


WDS 
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15.  Compute  the  output  power  of  DS 

NS 

P™  -  (NS)'  Z  P1k(0) 

CIS  J.k^i  JR 

and  signal -to -noise  ratio  gain  of  DS 


gds  =  10  logio  P, 


IN 


DS 


16.  Apply  the  filter  coefficients  to  obtain  the  FS  trace 
NFP-1 

2  NS 

s  -  y  y  a  x' 

NFP,n  NFp  !  V  Aj,m(NIM)+a,  n=  1,2,3,... 


17.  The  DS  and  WDS  traces  are  obtained  in  a  manner  similar  to  the  computation  in 
(16),  with  the  appropriate  weights  replacing  the  0jm's.  For  DS  the  weights  are 

0.  =  1/NS,  j  =  1,. . .  ,NS,  and  for  WDS  the  0.'s  are  as  computed  in  step  10. 

J  1 
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IX.  EXPERIMENTAL  RESULTS  AND  PRELIMINARY  EVALUATION  OF  FREQUENCY- 

DOMAEM  SYNTHESIS  PROCEDURE 


We  will  now  present  some  of  the  experimental  results  obtained  using  the 
frequency -domain  synthesis  computer  program.  Figures  1,  4,  7,  10,  13  show  the 
variation  of  anc*  signal -to -noise  ratio  gains  with  frequency.  The  original 

data  was  sampled  every  1/20  second  but  every  other  sample  was  used  in  the  fitting 
interval  so  that  the  sampling  rate  was  10  cps  corresponding  to  the  foldover  frequency 
of  5  cps  shown  in  the  figures.  Although  the  filters  are  designed  as  if  the  sampling 
interval  were  1/10  second,  the  filters  are  applied  to  the  data  with  the  sampling  interval 
of  1/10  second  but  in  steps  of  1/20  second.  In  Figs.  2,  5,  8,  11,  14  are  shown  the 
results  of  the  frequency-domain  processing.  The  top  trace  is  the  center  seismometer 
of  the  25-sensor  subarray,  the  next  trace  is  WDS,  the  next  is  DS,  and  the  bottom  trace 
represents  FS,  all  defined  previously.  The  corresponding  results  for  the  time -domain 
synthesis  are  shown  in  Figs.  3,  6,  9,  12,  15.  It  should  be  noted  that  only  the  FS 
trace  changes  between  these  two  sets  of  figures.  Similar  results  for  the  frequency- 
domain  synthesis  and  for  a  different  event  are  given  in  Figs.  16-24.  The  maimer  in 
which  the  signal -to -noise  ratio  gain  drops  outside  the  fitting  interval  for  the  time  and 
frequency -domain  methods  is  given  in  Fig.  25,  for  various  NFP  and  fitting  interval 
lengths.  The  fitting  interval  estimates  for  Pqut  were  compared  with  measured  values 
and  were  found  to  be  in  agreement  within  ±  1  db.  A  plot  of  the  frequency  window 
|W^(x)|  ,  and  —  10  log^  |  W^(x)/W^(0)  |  N=21,  used  in  the  spectral  estimation  program 

are  shown  in  Figs.  25  and  27,  respectively.  The  signal -to -noise  ratio  gain  vs.  frequency 
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was  obtained  with  somewhat  better  frequency  resolution  by  using  NFP  =  33,  as  indicated 
in  Fig.  28.  The  behavior  of  this  gain  at  low  frequencies  is  interesting  as  it  indicates 
that  the  gain  is  highest  at  about  0. 3  cps  and  drops  sharply  above  this  frequency. 

The  results  in  Figs.  16  through  24  for  the  2/1/66  Central  Kazakh  event  are 
rather  interesting  since  the  epicenter  for  this  event  is  located  about  900  miles  to  the 
southwest  of  a  presumed  nuclear  test  site.  The  raw  data  traces  do  not  show  clearly 
the  PcP  and  pP  phases,  as  can  be  seen  from  the  top  trace  in  these  figures.  However, 
the  FS  traces  do  show  these  phases  quite  clearly.  The  arrival  of  pP  at  about  seven 
seconds  after  the  initial  P  phase  establishes  the  focal  depth  of  the  event,  according  to 
seismological  travel  time  tables,  as  approximately  33  km.  This  in  turn  establishes  the 
event  as  being  definitely  an  earthquake.  This  conclusion  was  also  reached  by  the  U.S. 
Coast  and  Geodetic  Survey  on  the  basis  of  recordings  from  a  worldwide  network  of 
stations . 

The  experimental  results  show  that  the  frequency-domain  filter  tends  to  produce 
a  much  smaller  precursor,  in  general,  than  the  time-domain  filter.  This  can  be  seen 
by  comparing  the  results  on  the  11/11/65  Rat  Island  event.  In  addition,  the  signal -to- 
noise  ratio  gain  of  the  frequency-domain  filter  tends  to  remain  about  the  same  outside 
the  fitting  interval  as  it  is  inside  the  fitting  interval  for  periods  of  up  to  fourteen 
minutes  after  the  fitting  interval.  This  is  not  true  of  the  time -domain  filter,  which 
for  NFP  =  21  and  3  minute  fitting  intervals  tends  to  drop  about  4  db  just  outside  the 
fitting  interval.  However,  the  gain  of  the  time -domain  filter  is  still  about  2  db  better 
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than  that  of  the  frequency  domain  filter  outside  the  fitting  interval,  for  NFP  =21  and 
3  minute  fitting  interval.  This  loss  is  quite  acceptable  when  it  is  recognized  that  the 
7094  computer  running  time  to  synthesize  the  filter  is  only  about  4  minutes  for  the 
frequency -domain  method  as  opposed  to  30  minutes  for  the  time-domain  technique. 
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SJOB 
SIBSYS 
SATTACH 
SAS 

SATTACH 

SAS 

SEXECUTE 
SIBJOB 
SIBFTC  AA 

COMMON 


FREQ  DOMAIN  SYNTHESIS 


RJKOLKER 


B8 

SYSUT8 

A7 

SYSUT7 
IB  JOB 


/EDIT/ 
CALL  SPECTR 
CALL  FILTER 
CALL  DECOMP 
GO  TO  1 
END 


NYEAR*NDAY*NHR*NMIN*NSEC*MSC.NCH*ICH(255) 


SIBMAP 

UN21 

6 

ENTRY 

•UN21* 

•UN21* 

PZE 

UNIT21 

UNIT21 

FILE 

END 

F I LE2  *UT7  *  IN 

SIBMAP 

UN20 

6 

ENTRY 

•UN20* 

•UN20* 

PZE 

UNIT20 

UNIT20 

FILE 

END 

FI  LEI *UT8  *BI 

SIBMAP 

JUMP 

ENTRY 

•SKTBL 

•SKTBL 

PZE 

*+1. .LENGTH 

PZE 

•  UN20*  *  *20 

PZE 

•  UN21 *  *  *  21 

LENGTH 

EQU 

END 

♦-•SKTBL-1 

SORIGIN 

ALPHA 

SIBFTC 

SPCT 

‘256 


‘256 


C 

6002 


600 


500 


601 


SUBROUTINE  SPECTR 

DIMENSION  TBUF ( 255 ) *DBUF (132*25) 

COMPLEX  S ( 17 » 25 ) *CBUF (25*25) »F(17*325) *EXP  I (33*17) *C1 *C2 
COMMON  NS*NFP*NIM*AVE ( 25 ) *KD(25) *JCH(25) 

DIMENSION  IICH(25) »ISIES(25)*XC(25) * YC ( 25 ) 

COMMON/EDIT/  NYEAR *NDAY »NHR  *NMIN*NSEC*MSC*NCH* I CH( 255 ) 

EQUIVALENCE  ( CBUF (1*1) *DBUF ( 1 »1 ) ) 

INTEGER  WARMUP *BLKSZ*BUFSZ*FULLSZ.POCTR 

FORMAT ( 4H  NO* * 5X. 11HS I ESMOMETER * 5X . 1 1HCHANNEL  NO* » 5X . 1 1HCH*  ON  LIS 
IT  *  5X  *  1 1HXC00RDI NATE  *5X  * 11HYC00RDINATE  // ( 1 4,8X * A6 . 11X . I  5 *  1 IX  * 
2I5.5X.F11*3.5X.F11*3  >) 

WRITE! 6*600 ) 

FORMAT ( 45H1FREQUENCY  DOMAIN  SYNTHESIS  OF  FILTER  WEIGHTS  ) 

READ (5*500) NSKP .WARMUP .BLKSZtNSAMB  »NIM*NS*NFP  »POCTR 
FORMAT (81 5) 

NPRC* ( NSAMB*NIM*BLKSZ-1 ) /20 

WRITE(6»601) NSKP ♦ WARMUP  »BLXSZ *NSAMB*N IM  »NS*NFP  *NPRC*POCTR 
FORMAT (26H0NUMBER  OF 

*  44H  NUMBER  OF 

*  44H  NUMBER  OF 

*  44H  NUMBER  OF 

*  44H  NUMBER  OF 


RECORDS  SKIPPED  15// 

RECORDS  USED  GETTING  AVERAGES 
SAMPLE  POINTS  PER  SAMPLE  BLOCK 
SAMPLE  BLOCKS 
RECORDS  PER  SAMPLE  POINT 


15// 

15// 

15// 

15// 
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nr*  nn  n  on  nn  n  n  n  n  nn 


* 

* 


44H  NUMBER  OF  SIESMOMETERS 

46H  NUMBER  OF  FREOUENCI ES ( AND  FILTER  PTS/CHANNEL ) 
44H  NUMBER  OF  SECONDS  IN  FITTING  INTERVAL 
39H  FLAG  FOR  CONTROL  OF  MATRIX  PRINTOUTS 


15// 


15  // 
15  ) 


15// 


READ( 5*5001 M IICH< I ) • IS  I ES ( I ) *XC ( I) * YC ( I ) *  1*1 *NS) 

5001  FORMAT  C  I5*5X.A6*4X*2F10*3  ) 

REWIND  20 
REWIND  21 

READ( 20 )  NYEAR.NDAY*NHR*NMIN.NSEC.MSC*NCH*< ICH! I ) *1*1 *NCH) 

DO  1  1*1 *NS 

1  JCH< I ) *0 
IFLAG-O 

DO  2  I *1 *NS 
DO  2  J=1 *NCH 

IF(IICH(I).EQ.ICH(JJ )  JCH !  I  )  *  J 

2  CONTINUE 

DO  3  I *1 *NS 

IF ( JCH ( I ) *EO*  0 )  IFLAG*1 

IF (  JCH(I).EQ.  0)  WR I TE ( 6 » 6300 )  IICH(I) 

3  CONTINUE 

6300  FORMAT !  8H  CHANNEL  I4*23HN0T  ON  TAPE— -EXIT  TAKEN) 

IF ( I  FLAG*  ECU  1 )  STOP 

WRITE (6*6330)  NYEAR*NDAY»NHR»NMIN»NSEC»MSC 

6330  FORMAT!///  49H  DATING  INFORMATION  FOR  THE  EVENT  BEING  PROCESSED 
*//5H  YEARI4*3HDAYI4*4HHOURI4»3HMINI4*3HSECI4*3HMSCI4) 

WRITE! 6*6331)  NCH*( ICH! I ) *I»1*NCH) 

6331  FORMAT!  ////30H  INFORMATION  ABOUT  EDIT-3  TAPE  /// 

*30H  NUMBER  OF  CHANNELS  ON  TAPE  IS  I5.12HAND  THEY  ARE  / ! 20 1  5  )  ) 


WRITE (6 *6002) ( I  *  I  SI ES! I ) . I ICH 1 1 ) * JCH ( I ) *XC ! I ) *YC! I ) *  I  *  1 *NS ) 


NFPP* ! NFP-1 ) /2  +  1 

PI 2*2 *0*3 *141 5962 
YY«PI2/FLOAT!NFP-l) 


DO  4  L*1 *NFPP 

XL*L-1 

DO  4  M«1»BLKSZ 
FAC«XL*YY*FLOAT!M) 

4  EXPKM.L)  *  COS !  FAC  )*!1«0*0*0)+  S  IN !  FAC )  *  !0.0  *  1 .0 ) 


READ! 5*505 )  AZ*HVEL 
505  FORMAT ! 2F10*3 ) 

WRITE! 6*605 )  AZ»HVEL 

605  FORMAT!  7H  AZMUTH  F10*3*26H  HORIZONTAL  PHASE  VELOCITY  F10.3  ) 

CALL  DELAY!  AZ*HVEL*NS*XC*YC*KD) 

MI N* 100 
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non  nnn  on  on  on 


DO  5  I *1 tNS 

IF (  KD ( I ) «LT .MIN )  MIN*KD(I) 

5  CONTINUE 

DO  6  1*1. NS 

6  KD( I ) *KD( I ) -MIN 

WRITE! 6.606) ( KD( I ) .1*1 .NS) 

606  FORMAT ( 22H  BIASED  CHANNEL  DELAYS  // C  20 1  5  > I 


BUFSZ*NIM*BLKSZ 
FULLSZ*2*BUFSZ 
IF(FULLSZ.GT.  132)  IFLAG-1 
IFUFLAG.EQ.  1)  WRITE(6»607)  FULLSZ 

607  FORMAT (22H0D0UBLE  BUFFER  SIZE  OF  1 5 . 23HEXCEEDS  AVAILABLE  SPACE  ) 
I F ( IFLAG.EQ.l )  STOP 


CALL  AHEAD! 20. 0.NSKP) 

DO  7  1*1 *NS 

7  AVE! I ) *0 

DO  8  K*l. WARMUP 

READ! 20 ) ( TBUF ( I ) , I*1.NCH) 

DO  8  KS* 1 *NS 
N* JCH! KS ) 

8  AVE! KS ) *AVE ( KS ) +TBUF ( N ) 
XW*WARMUP 

DO  9  1*1. NS 

9  AVE! I ) *AVE ( I ) /XW 

WRITE! 6.702) (AVE! I ) .1*1. NS) 


DO  10  K-l.BUFSZ 
KK*BUFSZ+K 

READ (20) (TBUF! I ) .I*1.NCH) 

DO  10  KS*1.NS 
N* JCH ( KS ) 

10  DBUF ( KK .KS ) *TBUF ( N )-AVE ( KS ) 


DO  11  1*1.325 

DO  11  L* 1 .NFPP 
11  F(L.I)*(0. 0.0.0) 


DO  24  MAJ*1 .NSAMB 
DO  20  K-l.BUFSZ 
KK*BUFSZ+K 

READ (20) (TBUF! I ) » I*1.NCH) 

DO  20  KS'l.NS 

DBUF ( K .KS ) *DBUF ( KK .KS ) 

N* JCH ( KS ) 

20  DBUF ( KK .KS ) *TBUF ( N ) -AVE ( KS ) 
DO  22  LFR-l.NFPP 
DO  22  KS* 1 *NS 
C1*(0. 0.0.0) 

MT  *0 

DO  21  M*1 .BUFSZ .NIM 
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MT*MT+1 

MM*M+KD(KS) 

21  C1*C1+DBUF  (MM*KS)*EXPI  (MT *LFR  ) 

22  S(LFR*KS)=C1 

00  23  JS*1 *NS 

LSUM*JS*(JS-l>/2 
00  23  KS*1*JS 

LSM*LSUM+KS 
DO  23  LFR-l.NFPP 

23  F  ( LFR  *  LSM ) *S ( LFR  * JS )  *CON JG  ( S ( LFR  *KS )  )+F(LFR*LSM) 

24  CONTINUE 

XDI V-FLOAT  ( NSAMB )  *FLOAT  ( BLKSZ  ) 

00  25  1*1*325 

00  25  L*1*NFPP 

25  F(L*I )*F(L*I >/XDIV 
00  30  LFR*1 *NFPP 
DO  26  JS*1*NS 
LSUM*JS*(  JS-D/2 

DO  26  KS*1*JS 
LSM-LSUM+KS 

26  CBUF ( JS*KS ) *F ( LFR  * LSM ) 

00  27  JS*1.NS 

JJ*JS+1 

DO  27  KS* JJ»NS 

27  CBUF ( JS*XS ) *CONJG ( CBUF ( KS* JS )) 

WRITE (21 ) UCBUF(J*K) * J-l *NS ) *K*1 *NS ) 

IF  (  POCTR*NE.  0)  WRITE! 6*700 )  LFR 
IF (POCTR*EO*0 (  GO  TO  30 

00  28  K*1 *NS 

28  WR I TE (6  *701 )  ( CBUF ( J*K ) * J*1 *NS ) 

30  CONTINUE 

700  FORMAT ( 15H  PSD  MATRIX  FOR  I5*15H-1  TH  FREQUENCY  //) 

701  FORMAT (//( ( F12.4 *F12.4 ) * 7X ♦ ( F12 *4.F12.4 ) *7X* ( F12.4*F12 .4 ) . 7X. ( FI 1 
»*4*F12*4) ) ) 

END  FILE  21 
REWIN020 

WR I TE ( 21 )  ( AVE ( I ) *  1*1 *NS ) 

ENO  FILE  21 

702  FORMAT ( 47H  AVERAGE  DC  LEVELS  PREVIOUS  TO  FITTING  INTERVAL  //(12 
*F10.4) ) 

REWIND  21 

RETURN 

ENO 

SIBFTC  OLAY 

SUBROUTINE  DELAY ( AZ*HVEL*NS*XC*YC.KO) 

DIMENSION  XC(NS) *YC(NS) *KD(NS) *RS( 100) *TH( 100) 

00  3  1*1 *NS 

IF (SORT (XC( I )**2+YC( I ) **2 ) *EQ*0*0 )  GO  TO  2 
RSU)*SQRT(XCU  )**2  +YC(I)#*2) 

TH(I)*ATAN2(  YC(I)*XC(I)) 

GO  TO  3 

2  TH( I ) *0*0 
RS( I )*0*0 

3  CONTINUE 
C0NST*180. 0/3*14159 
DO  66  1*1 *NS 

66  TH( I )*TH( I )*CONST 
T*0*05 

C*3* 1415927/180*0 
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CC*-1*0/(HVEL*T) 

DO  50  I *1 #NS 

DEL  *RS( I)*CC*COS(C*(TH( I )-AZ) ) 

IF (  DEL.LE.0.0 )G0T025 
60  TO  30 

25  KD< I )*DEL-0*5 
GO  TO  50 

30  KD(I)*DEL+0#5 
50  CONTINUE 
RETURN 
END 

$0R IGI N  ALPHA 

SIBFTC  SOSO  FULIST 

SUBROUTINE  FILTER 

COMPLEX  CBUF(25*25) t  ABUF ( 25 ) *QBUF (25*25) *C1*C2*C3 
DIMENSION  W<25*33) ,OPF(25) #RHO( 26*26) *RR( 26.26) 

COMMON  NS *NFP*NIM * AVE ( 25 ) *KD ( 25 ) . JCH ( 25 ) 

6800  FORMAT (  21H  VARIANCE  OF  INPUT  IS  F13.4.5X* 

121HVARIANCE  OF  OUTPUT  IS  F13.4/ 

225H  NOISE  REDUCTION  FOR  WS.  F13*4*5X* 

♦24HNOISE  REDUCTION  FOR  WDS  F13.4  // 

*20H  VARIANCE  OF  WDS  F13.4  ) 

6900  FORMAT ( 4H  FORI  5  .25HTH  MULTIPLE  OF  FREQUENCY  /  19H  INPUT  VARIANC 
IE  IS  F13*4* 18HOUTPUT  VARIANCE  IS  F13.4.25HNOISE  REDUCTION  IN  DB  I 
2S  F13*4  ) 

6802  FORMAT { 23H  OPTIMAL  ONE  PT  FILTER  //(9F13.4)) 

6801  FORMAT ( 13H  FOR  CHANNEL  I5»18HFILTER  WEIGHTS  ARE// 

1  (9F14.8)) 

SUMO-O 

NFPP* ( NFP— 1 ) /2  +  1 
XNS*NS 
NFPM*NFP-1 
XNF*NFPM 

XDEN* I XNF+1 *0 ) / ( XNF*XNF ) 

REWIND  21 

TRIG*  2.0*3*141592/XNF 
SUMI *0 

XDI V*XNS*XNF 

DO  30  LFREQ*1 *NFPP 

LF*LFREQ-1 

READ (21 ) ( ( CBUF ( I t J) > 1*1 *NS) • J«1 >NS) 

XLF-LF 
Cl-<0. 0,0.0) 

DO  10  J* 1 *NS 
10  C1*C1+CBUF ( J  * J ) 

VI*REAL(C1)/XDIV 
C  NOW  THE  COMPLEX  INVERSION 

CALL  COMINV  ( 25*NS.CBUF.QBUF ) 

C1*(0*0*0.0) 

DO  14  J* 1 *NS 
DO  14  I *1 »NS 

14  C1*C1+0BUF( I *J) 

VO*REAL (Cl ) 

VO*XDEN/VO 

DO  16  K*1 *NS 

C2*(0,0*0.0) 

DO  15  J*1 *NS 

15  C2*C2+QBUF(K.J) 

16  ABUF  ( 1C)  *C2/C1 
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IF ( LFREQ.EQ*  1 )  60  TO  18 
IF(LFREQ.EQ.NFPP)GO  TO  21 
DO  17  J*1 »NS 

FAC1*REAL ( ABUF ( J ) ) 

FAC2*AIMAG( ABUF (J)) 

DO  17  K«1.NFP 

XXK=K-1 

YYK»-FLOAT(  (NFP-1  )/2M-XXK 

TAC»FAC1*C0S(YYK*XLF*TRIG)+FAC2*SIN!YYK*XLF*TRIG) 

17  W( J.K)*W( J.K)+2.0*TAC 
GO  TO  23 

18  DO  19  J«1.NS 

DO  19  K*1 *NS 

19  RHO! K* J ) =REAL ( CBUF ( X* J ) ) 

DO  20  J=1 *NS 

FAC1«  REAL(ABUF( J) ) 

DO  20  K-l.NFP 

20  W(J.M«FAC1 
GO  TO  23 

21  DO  22  J=1 *NS 

DO  22  K=1 *NFP 

KK*K-1 

KZ«KK-(NFP-l)/2 
KZ* I ABS! KZ ) 

FAC1«1.0 

DO  2022  11=1. NFP 
IF(II-l-KZ)  2022.22.2022 
2022  FAC1—FAC1 

22  W( J.K)«W(J.K)+FAC1*REAL(ABUF(J) ) 

DO  222  J*1.NS 

DO  222  K*1.NS 

222  RHO(K.J)*RHO!K,J)+  REAL  (CBUF(K.J)) 

23  IF(LFREQ.EQ.l)  GO  TO  25 
IF(LFREQ.EQ.NFPP)  GO  TO  25 
DO  24  J*1.NS 

DO  24  K*1.NS 

24  RHO( J»K)=RHO( JtK)t  2.0*REAL ( CBUF ( J .K )) 

25  IF  (LFREQ.EQ. 1)  GO  TO  26 

IF ( LFREO.EQ.  NFPP )  GO  TO  26 
SUMO*SUMO+2.0»VO 
SUMI * SUM 1+2 .0*V  I 
GO  TO  28 

26  SUMO-SUMO+VO 
SUMI *SUMI+VI 

28  DBR*  10.0*AL0610(VI/V0) 

WRITE! 6.6900 )  LF. VI. VO. DBR 

30  CONTINUE 

CALL  AHEAD (21.2.0) 

DO  31  K*1.NFP 

DO  31  J«1.NS 

31  W( J»K)*W(J.K)/FLOAT(NFP-l) 

WRITE! 21 ) ( <W( J.K) *J*1»NS> .K*l .NFP) 

DO  32  J-l.NS 

DO  32  K-l.NS 

32  RHO! X » J ) *RHO( K. J ) /FLOAT (NFP ) 

NS1*NS+1 

DO  33  I =1 .NS 
RHO! NS1 *  I ) *1 .0 

33  RHO! I *NS1 )=1.0 
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n  * 


RHO( NS1 *NS1 ) “0 

CALL  MAT INV! 26.NS1 .RHO.RR ) 

DO  34  I-l.NS 

34  OPF( I ) -RR ( I *NS1 ) 

WRITE(21 ) (OPF ( I ) .I-l.NS) 

ENDFILE  21 

REWIND  21 

DBR- 10 • 0*AL0G1 0 ( SUM I /SUMO ) 

DBS*10*0*ALOG10(-SUMI /RR C  NS1 *NS1 ) ) 

VAR--RR ( NS1 *NS1 ) 

WRITE! 6*6800 )  SUMI .SUMO. DBR .DBS .VAR 
WRITE  C  6.6902 ) 

6902  FORMAT  C 36H  ZERO  SHIFT  CROSS-CORRELATION  MATRIX  ) 

DO  333  I-l.NS 

333  WRITE! 6.6901 )  (RHO! I . J) . J-l.NS) 

6901  FORMAT (//(9F13.4) ) 

DO  35  J-l.NS 

35  WRITE! 6.6801 )  J, (W( J.K) .K-l.NFP) 

WRITE! 6.6802 )  I OPF ! I ) . I-l.NS ) 

DO  356  K-l.NFP 
SUM-0 

DO  355  J-l.NS 

355  SUM-SUM+W ( J  »  K  ) 

356  WRITE (6.6922 )  K.SUM 

6922  FORMAT ( 29H  SUM  OF  WEIGHTS  AT  FILTER  PT#  I3.7X.  F20.9) 
RETURN 
END 

SIBFTC  CINV 

SUBROUTINE  COMINVtNDIM.NORD.C.CI ) 

COMPLEX  C ( NDIM.NDIM ) .Cl (NDIM.NDIM) 

DIMENSION  BI50.50) »BI (50.50) 

DO  10  J-l.NORD 
JJ-J+NORD 
DO  10  I-l.NORD 
I I-I+NORD 

B ( I . J ) -REAL ( C ( I » J ) ) 

B ( I . JJ ) -A  I MAG (C ( I ♦ J ) ) 

B( 1 1 . JJ ) -REAL (C ( I *J) ) 

B ( 1 1 . J ) --AIMAG ( C ( I » J ) ) 

10  CONTINUE 

NDORD«2*NORD 

CALL  MATINV! 50.NDORD.B.BI ) 

DO  20  J-l.NORD 
JJ* J+NORD 
DO  20  I-l.NORD 

20  CI(I»J)«BI(I.J)*(1.0.0.0)+BI(I.JJ)*(0.0»1.0) 

RETURN 

END 

IBFTC  MAIN 

MATRIX  INVERSION  BY  BORDERING 
SUBROUTINE  MATINV(NU*NORD.A.AI ) 

DIMENSION  A(NU.NU) *AI (NU.NU) 

DIMENSION  I CH ( 50 ) *KCH ( 50 ) .COL! 50) .ROW! 50 ) 

DO  1  K-l.NORD 
(CCH(K)-X 

1  ICH! K) -X 

DO  2  K-l.NORD 
IF  ( A(K.l ) )  3.2,3 

2  CONTINUE 


46 


WRITE( 6*1000 ) 

1000  FORMAT ( 1H  *22H  EXIT1*MATRIX  SINGULAR) 
CALL  EXIT 

3  ICH( 1 ) *X 

DO  4  J-l*NORD 
TEMP-A(1*J) 

A< 1#J)«A<X.J) 

4  A( X* J) -TEMP 

AI ( 1 *1 )  •  1*0/A (1*1) 

DO  15  N-2  *NORD 
DO  111  NX*N#NORD 
XCH{ N ) -NX 
DO  90  JB1 *NORD 
TEMP»A(J*NX) 

A ( J  t  NX ) - A  C  J  #  N ) 

90  A( J«N) -TEMP 

DO  11  XX-N»NORD 
ICH(N)«XX 
DO  5  J*1 »NORD 

TEMP-A (  N  *  J  ) 

A ( N* J ) -A ( XX ♦ J ) 

5  A  <  XX  t J )  “TEMP 
L-N-l 

DO  7  I *1  *L 
SUM1-0.0 
SUM2-0.0 
DO  6  J«1*L 

SUM1-SUM1+AI ( I#J)#A(J,N) 

6  SUM2*SUM2+A( N* J )*AI ( J*  I ) 

COL ( I ) -SUM1 

7  ROW( I ) -SUM2 
SUM1-0.0 
SUM2-0.0 

DO  8  J-l  »L 

SUM1-  SUM1  +  A(N*J)*COL< J) 

8  SUM2-  SUM2+  A ( J *N ) *ROW ( J ) 

AVE  *  ( SUM1  +  SUM2 ) /2«0 
ADIV-  A( N*N)  -AVE 

IF (ADI V)  12  *9#  12 

9  ICH( N) *N 

DO  10  J*1 tNORD 
TEMP-  A ( N* J ) 

A ( N* J ) -A ( XX  * J ) 

10  A ( XX ♦ J ) -TEMP 

11  CONTINUE 
XCH(N)*N 

DO  91  J*1 »NORD 
TEMP-A ( J»NX) 

A( J*NX ) -A ( J*N ) 

91  A ( J*N ) -TEMP 
111  CONTINUE 

WRI TE ( 6  » 1001 ) N 

1001  FORMAT ( 1H  *21HEX I T2 * MATR IX  SINGULAR 
CALL  EXIT 

12  AI(N*N)-1*0/ADIV 
DO  13  J*1 *L 

AI  ( J»N)— COL<  JJ/ADIV 

13  AI  <N#J)—ROW<  JJ/ADIV 
DO  14  I-ltL 
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DO  14  J  =  1  *L 

14  AI ( I *J)*AI ( I*J)+COL( I )*ROW( JJ/ADIV 

15  CONTINUE 

DO  116JJ*1*N0RD 
KR*NORD+l-JJ 
KL«KCH ( KR ) 

DO  116  J=1 *NORD 
TEMP»A(J*KR) 

A( J*KR)»A( J.KL) 

A ( J.KL )  “TEMP 
TEMP-AI <KR*J) 

A I  ( KR  »  J ) * A I ( KL  »  J  ) 

116  AI (KL*J)*TEMP 
DO  16  J*1 *NORD 
KR=NORD+l-J 
KN*ICH (KR ) 

DO  16  LL=1 *NORD 
TEMP*A(KR.LL) 

TEMPP»AI ( LL  *KR ) 

A(KR.LL)  *A< KNtLL ) 

AI  (LL.KR)*AI (LL*KN) 

A(KN*LL)=TEMP 

16  AI (LL*KN)«TEMPP 
RETURN 

END 

SORIGIN  ALPHA 

SIBFTC  PROC  LIST 

SUBROUTINE  DECOMP 

COMMON  NS*NFP»NIM»AVE ( 25 ) »KD ( 25 ) » JCH ( 25 ) 

COMMON  /EDIT/  NYEAR*NDAY »NHR»NMIN#NSEC*MSC»NCH* ICH( 255 ) 
COMMON/BETH/  NJMP*NSQSM  »W ( 25 *33 ) *OPF ( 25 ) 

READ ( 5*5000 )  NJMP *NSQSM*NSKP .NPLOT 
5000  FORMAT (415) 

REWIND  20 
REWIND  21 

READ( 20 )  NY*ND*NH.NMN.NSC*MSC*NCH*( ICH( I ) *I=1*NCH) 

CALL  AHEAD ( 21*2*0) 

READ ( 21 M (W( J*K) *J=1.NS) .K-l.NFP) 

READ (21) (OPF(K) *K*1*NS) 

IF(NJMP*EO.O)  GO  TO  10 
CALL  PLOTT 

10  IF(NSKP.EQ.O)  GO  TO  20 
20  RETURN 
END 


SORIGIN 

BETA 

$  JOB 

SIBSYS 

SEXECUTE 

IB  JOB 

SIBJOB 

SIBFTC  PLTT 

FUL 1ST 

SUBROUTINE  PLOTT 

COMMON  NS.NFP *NIM * AVE (25) *KD(25) » JCH ( 25 ) 

COMMON  /EDIT/  NYEAR *NDAY *NHR ,NMI N. NSEC. MSEC.NCH • ICH ( 255 ) 
COMMON  /BETH/  NJMP*NSOSM*W( 25*33) »OPF( 25 ) 

DIMENSION  TBUF ( 255 ) *DBUF (25*200) *IYLINE(4) *XSEC ( 6 ) » I XSEC ( 6 ) 
DIMENSION  MB ( 4 ) *MT ( 4 ) 

DIMENSION  ARRAY( 1000*4) 

DIMENSION  LD ( 25 ) 

INTEGER  CHOP. DELTA 
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NDEL= { NFP— 1 ) *NIM/2 
DO  970  1*1. NS 
LD ( I ) »KD( I ) 

970  KD( I ) *KD( I )+NDEL 

WRITE(6*8700)  NJMP.NSQSM 
8700  FORMAT (26H  NUMBER  OF  RECORDS  SKIPPED  I5t 
*  29H  NUMBER  OF  RECORDS  PROCCESSED  15  ) 

READ (5.570)  TRCMX.NSITE 
570  FORMAT (  F10.3,  15) 

TRCMN*-TRCMX 

CALL  STOIDV! 15HFOUR  TRACE  PLOT  .3) 

CALL  BRITEV 
CALL  CAMRAV! 9 ) 

CALL  FRAMEVCO) 

CALL  PRINTV!-13.13HDATE  OF  EVENT  .400.670) 

CALL  PRINTV  (-4.4HYEAR  .450.630) 

CALL  PRINTV(-3.  3HDAY  .450*600) 

CALL  PRINTV ( -21 .  21HSTARTING  TIME  OF  PLOT  .400.560) 
CALL  PRINTV  (  -4  .4HH0UR  .450.530) 

CALL  PRINTV ( -6  .6HMINUTE  .450.500) 

CALL  PRINTV! “6. 6HSEC0ND  .450.470) 

CALL  PRINTV (-11.11HSITE  NUMBER  .400.420) 

CALL  PRINTVI-2.2HFS  .950.148  ) 

CALL  PRINTV1-3.3HWDS  .950.648) 

CALL  PRINTV ! -2.2HDS  .950.398) 

CALL  PRINTVI-4.4HDEEP  .950.898) 

CALL  PRINTV C -4.4HWELL  .950.882) 

DELTA»50*NSQSM 

CH0P«MSC+1000*NSEC+60000*NMIN+3600000*NHR 
CH0P=DELTA+CH0P 
I HR*CHOP/ 3600000 

CHOP* CHOP-3600000* ( CHOP/3600000 ) 

IMIN*CHOP/60000 
CHOP«CHOP-60000* ! CHOP/60000 ) 

ISEC-CH0P/1000 

X-NYEAR 

CALL  LABLV (X. 530. 630. 4. 1.4) 

X*NDAY 

CALL  LABLVIX.530.600.4.1.4) 

X*  IHR 

CALL  LABLV (X.530.530.4.1.4) 

X*IMIN 

CALL  LABLV ( X. 530. 500 .4. 1 .4 ) 

X*ISEC 

CALL  LABLV1X.530.470.4.1.4) 

X*NSITE 

CALL  LABLV! X. 530.420 .4. 1 .4 ) 

CALL  FRAMEVIO) 

I YLINE 1 1 ) *148 
I YLI NE ( 2 ) *648 
I YL INE ( 3 ) *398 
I YLINE !4)»898 
XSEC ( 1 ) *0 
XSEC ( 2 ) *10.0 
XSEC ( 3 ) *20.0 
XSEC ( 4 ) *30.0 
XSEC ( 5 ) *40.0 
XSEC ! 6 ) *50.0 
IXSEC ( 1 ) *12 
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IXSEC(2>*212 

IXSECt 3)*412 

IXSEC ( 4 ) =612 

IXSEC ( 5 ) *812 

I XSEC ( 6 ) *978 

00  30  1*1*4 

MT( I )*1023-( IYLINE( I) +125) 

30  MB( I )*IYLINE( I )-125 
CFS*0 
CWDS*0 
CAP-0 
XNS-NS 
KBL-0 

WRI TE ( 6*6300 ) 

WRI TE( 6*6400 ) 

NN-NJMP+1 

CALL  AHEAD  (20*0«NN) 

DO  35  1*101*200 

READ (20) (TBUF(K) *K«1*NCH> 

DO  35  KS*1 *NS 
N= JCH ( KS ) 

35  DBUF  ( KS*  I )  *TBUF  ( N  )-AVE  ( KS ) 
NBLOK*NSQSM/100 

50  IF(NBLOK*LE*0)  GO  TO  71 

IF ( N8LOK-10 )  51*52*52 

51  IBMAX-NBLOK 
GO  TO  53 

52  IBMAX-10 

53  NBLOK-NBLOK-IO 

DO  58  IB*1 • IBMAX 
KBL-KBL+1 
DO  54  1*1*100 

11*1+100 

READ (20) (TBUF(X) »K*1*NCH) 

DO  54  KS=1 *NS 
DBUF ( KS  *  I ) *DBUF ( KS • 1 1  ) 

N-JCH(KS) 

54  DBUF ( KS* 1 1 ) *TBUF ( N ) -AVE ( KS ) 

I ADX* ( IB-1 ) *100 

BFS-0 

BWDS-0 

BAP-0 

DO  57  I  T*1 *  100 

SUMAP-0 

SUMC*0 

SUMW-0 

SUMD-0 

IADX-I ADX+1 

DO  56  ICS-1  *NS 

IDFS* I T+LD( KS ) -NIM 

IDWS* IT+XD( XS ) 

SUHF-0 

DO  55  KF-l.NFP 
IDFS* IDFS+NIM 

55  SUMF-SUMF+DBUF  ( KS  •  IDFS  )*W  ( KS*ICF  ) 
SUHC-SUMC+SUMF 

SUMW«SUMW+DBUF ( KS  *  I DWS ) *OPF ( KS ) 
SUMD-SUMD+DBUF ( KS  » I DWS ) 

56  SUMAP-SUMAP  +  ( DBUF ( KS  *  I DWS ) **2 ) 


50 


BFS=BFS+SUMC*SUMC 
BWDS=BWDS+SUMW*SUMW 
BAP=BAP+ ( SUMAP /XNS ) 

ARRAY ( I ADX • 1 ) *$UMC 
ARRAY ( I ADX .2 ) *SUMW 
ARRAY ( I ADX  *  3 ) *SUMD/XNS 
I  Z=  I T+KD ( 1 ) 

57  ARRAY< IADX*4)*DBUF(1.IZ) 

BFS=BFS/100.0 

BWDS=BWDS/100.0 

BAP*BAP/100.0 

CFS*CFS+BFS 

CAP*CAP+BAP 

CWDS*CWDS+BWDS 

58  WR I TE < 6*6200 ) KBL * BFS ♦ CFS * BWDS *CWDS .BAP .CAP 
I AD* I ADX-1 

1X1*12 
1X2*1012 
DO  60  1*1.4 

60  CALL  L1NEV( IX1»IYLINE( I ) .IX2.IYLINEI I ) ) 

DO  61  1*1.6 

61  CALL  LABLV(XSEC<n.IXSEC(  I)  .12.4.1.4) 

DO  65  1*1.6 

65  XSECU  )*XSEC<  I  M-50.0 
DO  62  1*1.4 

IXS* 12 

CALL  YSCALVi TRCMN*TRCMX.MB( I ) .MT ( I ) ) 

DO  62  KP*1 *  I AD 
IXL* IXS+1 
Z1*ARRAY(KP.I ) 

Z2* ARRAY ( KP+1 » I ) 

IY1*IYV(Z1) 

I Y2* I YV ( Z2 ) 

CALL  LINEV( IXS.IY1.IXL.IY2) 

62  IXS* IXL 

CALL  FRAMEV(O) 

60  TO  50 
71  CALL  PLTND 
REWIND  20 
RETURN 

6400  FORMAT ( 7H  BLOCK  * 12X.2HFS. 14X .7HFS  CUM. .14X.3HWDS. 13X .8HWDS  CUM.. 

S15X.2HAP .14X.7HAP  CUM.  //) 

6200  FORMAT (I6.6(7X.F12.4)) 

6300  FORMAT (14H  SQASM  OUTPUT  //////) 

END 
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3-64-4959 


Figure  1.  Input  power,  output  power  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequency -domain  maximum -likelihood  filter, 
11/11/65  Rat  Island  event,  subarray  Bl,  NIM  =  2,  NFP  =  21, 
3  minute  fitting  interval. 
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3-64-4960 


Figure  2.  Processed  traces  obtained  by  frequency -domain  synthesis 

procedure,  11/11/65  Rat  Island  event,  subarray  Bl,  NIM  =  2, 
NFP  =  21,  3  minute  fitting  interval. 
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Figure  3.  Processed  traces  obtained  by  time-domain  synthesis  procedure, 

11/11/65  Rat  Island  event,  subarray  Bl,  NIM  =  2,  NFP  =  21, 

3  minute  fitting  interval. 


3-64-4962 


Figure  4.  Input  power,  output  power  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequency-domain  maximum -likelihood  filter, 
11/11/65  Rat  Island  event,  subarray  A(J,  NIM  =  2,  NFP  =  21, 
3  minute  fitting  interval. 


56 


3-64-4963 
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Figure  5.  Processed  traces  obtained  by  frequency-domain  synthesis 

procedure,  11/11/65  Rat  Island  event,  subarray  A0,  NIM  =  2, 
NFP  =  21,  3  minute  fitting  interval. 
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Figure  6. 


Processed  traces  obtained  by  time-domain  synthesis  procedure, 
11/11/65  Rat  Island  event,  subarray  A0,  NIM  =  2,  NFP  =  21, 

3  minute  fitting  interval. 
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Figure  7.  Input  power,  output  power  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequenoy -domain  maximum -likelihood  filter, 
11/11/65  Rat  Island  event,  subarray  B3,  NIM  =  2,  NFP  =  21, 
3  minute  fitting  interval. 
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Figure  8.  Processed  traces  obtained  by  frequency-domain  synthesis 

procedure,  11/11/65  Rat  Island  event,  subarray  B3,  NIM  =  2, 
NFP  =  21,  3  minute  fitting  interval. 
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Figure  9.  Processed  traces  obtained  by  time-domain  synthesis  procedure, 

11/11/65  Rat  Island  event,  subarray  B3,  NIM  =  2,  NFP  =  21, 

3  minute  fitting  interval. 
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Figure  10.  Input  power,  output  power  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequency -domain  maximum -likelihood  filter, 
11/11/65  Rat  Island  event,  subarray  B4,  NIM  =  2,  NFP  =  21, 
3  minute  fitting  interval. 
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Figure  11.  Processed  traces  obtained  by  frequency-domain  synthesis 

procedure,  11/11/65  Rat  Island  event,  subarray  B 4,  NIM  =  2, 
NFP  =  21,  3  minute  fitting  interval. 
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Figure  12.  Processed  traces  obtained  by  time -domain  synthesis  procedure, 
11/11/65  Rat  Island  event,  subarray  B4,  NIM  =  2,  NFP  =21, 

3  minute  fitting  interval. 
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Figure  13.  Input  power,  output  power  and  signal -to-noise  ratio  gain  vs. 

frequency  for  frequency-domain  maximum -likelihood  filter, 
11/11/65  Rat  Island  event,  subarray  B2,  N1M  =  2,  NFP  =  21, 
3  minute  fitting  interval. 
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Figure  14.  Processed  traces  obtained  by  frequency-domain  synthesis 

procedure,  11/11/65  Rat  Island  event,  subarray  B2,  NIM  =  2, 
NFP  =  21,  3  minute  fitting  interval. 
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Figure  15.  Processed  traces  obtained  by  time -domain  synthesis  procedure, 
11/11/65  Rat  Island  event,  subarray  B2,  NIM  =  2,  NFP  =  21, 

3  minute  fitting  interval. 
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Figure  16.  Input  power,  output  power  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequency-domain  maximum -likelihood  filter, 
2/1/66  Central  Kazakh  event,  subarray  Bl,  NIM  =  2,  NFP  =  21, 
3  minute  fitting  interval. 
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Figure  17.  Processed  traces  obtained  by  frequency-domain  synthesis 
procedure,  2/1/66  Central  Kazakh  event,  subarray  Bl, 
NIM  =  2,  NFP  =  21,  3  minute  fitting  interval. 
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Figure  18.  Input  power,  output  power  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequency-domain  maximum -likelihood  filter, 
2/1/66  Central  Kazakh  event,  subarray  A0,  NIM  =  2, 

NFP  =  21,  3  minute  fitting  interval. 
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Figure  19.  Processed  traces  obtained  by  frequency -domain  synthesis 
procedure,  2/1/66  Central  Kazakh  event,  subarray  A0, 
NLM  =  2,  NFP  =  21,  3  minute  fitting  interval. 
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Figure  20.  Input  power,  output  power,  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequency-domain  maximum -likelihood  filter, 
2/1/66  Central  Kazakh  event,  subarray  B3,  NIM  =  2,  NFP  =21, 
3  minute  fitting  interval. 
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Figure  21.  Processed  traces  obtained  by  frequency-domain  synthesis 
procedure,  2/1/66  Central  Kazakh  event,  subarray  B3, 
NIM  =  2,  NFP  =  21,  3  minute  fitting  interval. 
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Figure  22.  Input  power,  output  power  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequency-domain  maximum  -likelihood  filter, 
2/1/66  Central  Kazakh  event,  subarray  B4,  NIM  =  2,  NFP  =  21, 
3  minute  fitting  interval. 
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Figure  23.  Processed  traces  obtained  by  frequency-domain  synthesis 
procedure,  2/1/66  Central  Kazakh  event,  subarray  B4, 
NIM  =  2,  NFP  =  21,  3  minute  fitting  interval. 
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Figure  24.  Processed  traces  obtained  by  frequency-domain  synthesis 
procedure,  2/1/66  Central  Kazakh  event,  subarray  B2, 
NIM  =  2,  NFP  =  21,  3  minute  fitting  interval. 
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Figure  25.  Signal -to -noise  ratio  gain  of  maximum -likelihood  filter 

outside  the  fitting  interval  for  various  NFP  and  fitting 
interval  lengths,  noise  from  2/4/66,  subarray  AfJ, 
Azimuth  =  0°,  Horizontal  Velocity  =  15  km /sec,  SNR 
gain  measured  over  2  minute  intervals. 
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Figure  26.  Plot  of  frequency  window  [W^(rrf/5)  ]  ,  N  =  21. 
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Figure  27.  Hot  of  frequency  window  -  10  log._  [  W  (TTf/5)/W  (0)]2, 
N  =  21.  1U  W  N 
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Figure  28.  Input  power,  output  power  and  signal -to -noise  ratio  gain  vs. 

frequency  for  frequency -domain  maximum -likelihood  filter, 
11/11/65  Rat  Island  event,  subarray  Bl,  NIM  =  2,  NFP  =  33, 
3  minute  fitting  interval. 
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